Return to Main Page

\(\color{darkblue}{\textbf{sf}}\)


Used for processing spatial data and basic visualization

library(tidyverse)
library(sf)
library(raster)

\(\color{dodgerblue}{\textbf{Importing}}\)

\(\color{skyblue}{\textrm{- Reading in Data}}\)

Vector data (in any spatial form) can be read in with read_sf…

url <- "https://github.com/social-lorax/new_guides/raw/main/Spatial/Borough%20Boundaries.geojson"
boros <- read_sf(url) 


…or transformed using st_as_sf()

spatial <- read_csv("pathway/to/spatial.csv") %>% 
  st_as_sf()


Raster data is read in using raster() if it is single-band…

url <- "https://github.com/social-lorax/new_guides/blob/main/Spatial/canopy.tif?raw=true"
canopy <- raster(url)


…and brick() if it is multi-band

url <- "https://github.com/social-lorax/new_guides/blob/main/Spatial/manhattan.tif?raw=true"
manhattan <- brick(url)

\(\color{dodgerblue}{\textbf{Removing Geometry}}\)

Having a geometry can be annoying for other functions, so st_set_geometry(NULL) removes it

non_spatial <- spatial %>% 
  st_set_geometry(NULL)

\(\color{dodgerblue}{\textbf{Geocoding}}\)

boro_halls <- tibble(hall = c("New York City Hall", "Bronx County Courthouse", "Brooklyn Borough Hall",
                              "Queens Borough Hall", "Staten Island Borough Hall"),
                     street = c("New York City Hall, New York, NY 10007", 
                                "851 Grand Concourse, Bronx, NY 10451", 
                                "209 Joralemon St, Brooklyn, NY 11201",
                                "120-55 Queens Blvd, Queens, NY 11424", 
                                "10 Richmond Terrace, Staten Island, NY 10301"))

boro_halls
hall street
New York City Hall New York City Hall, New York, NY 10007
Bronx County Courthouse 851 Grand Concourse, Bronx, NY 10451
Brooklyn Borough Hall 209 Joralemon St, Brooklyn, NY 11201
Queens Borough Hall 120-55 Queens Blvd, Queens, NY 11424
Staten Island Borough Hall 10 Richmond Terrace, Staten Island, NY 10301


\(\color{skyblue}{\textrm{- Cleaning Messy Addresses}}\)

library(postmastr) #remotes::install_github("slu-openGIS/postmastr")

#Set-up
boro_halls <- boro_halls %>% 
  pm_identify(var = "street")

boro_halls_unique <- boro_halls %>% 
  pm_prep(var = "street", type = "street")


ZIP Code

#Check
boro_halls_unique %>% pm_postal_none()
pm.uid pm.address
#Parse
boro_halls_unique <- boro_halls_unique %>% 
  pm_postal_parse()


State

#Check
nyDict <- pm_dictionary(locale = "us", type = "state", filter = "NY", case = c("title", "upper"))

boro_halls_unique %>% pm_state_none(dictionary = nyDict)
pm.uid pm.address pm.zip
#Parse
boro_halls_unique <- boro_halls_unique %>% 
  pm_state_parse()


City

#Check
defaultDict <- pm_dictionary(locale = "us", type = "city", filter = "NY")

specificDict <- pm_append(type = "city",
                          input = c("New York", "New York City", "NYC"     , "Bronx", "The Bronx",
                                    "Brooklyn", "Queens", "Jamaica", "Staten Island"),
                          output = c(NA       , "New York"     , "New York", NA     , "Bronx"    , 
                                     NA       , NA      , "Queens" , NA))

boro_halls_unique %>% pm_city_none(dictionary = specificDict)
pm.uid pm.address pm.state pm.zip
#Parse
boro_halls_unique <- boro_halls_unique %>% 
  pm_city_parse(dictionary = specificDict)


House Number

#Check
boro_halls_unique %>% pm_house_none()
pm.uid pm.address pm.city pm.state pm.zip
1 New York City Hall New York NY 10007
#Parse after fixing or ignoring
boro_halls_unique <- boro_halls_unique %>% 
  pm_house_parse()


Street Direction (e.g., East Broadway)

#Check 
boro_halls_unique %>% pm_streetDir_none()
pm.uid pm.address pm.house pm.city pm.state pm.zip
1 New York City Hall NA New York NY 10007
2 Grand Concourse 851 Bronx NY 10451
3 Joralemon St 209 Brooklyn NY 11201
4 Queens Blvd 120-55 Queens NY 11424
5 Richmond Terrace 10 Staten Island NY 10301

Moving on since none have a prefix


Street Suffix (e.g., St., Ave.)

#Check
boro_halls_unique %>% pm_streetSuf_none()
pm.uid pm.address pm.house pm.city pm.state pm.zip
1 New York City Hall NA New York NY 10007
2 Grand Concourse 851 Bronx NY 10451
#Parse after fixing or ignoring
boro_halls_unique <- boro_halls_unique %>% 
  pm_streetSuf_parse()


Street

#Street is assumed 
boro_halls_unique <- boro_halls_unique %>% 
  pm_street_parse()


Reassemble

boro_halls_parsed <- pm_replace(boro_halls_unique, source = boro_halls) %>% 
  pm_rebuild(output = "full", keep_parsed = "no")

boro_halls_parsed
hall street pm.address
New York City Hall New York City Hall, New York, NY 10007 New York City Hall New York NY 10007
Bronx County Courthouse 851 Grand Concourse, Bronx, NY 10451 851 Grand Concourse Bronx NY 10451
Brooklyn Borough Hall 209 Joralemon St, Brooklyn, NY 11201 209 Joralemon St Brooklyn NY 11201
Queens Borough Hall 120-55 Queens Blvd, Queens, NY 11424 120-55 Queens Blvd Queens NY 11424
Staten Island Borough Hall 10 Richmond Terrace, Staten Island, NY 10301 10 Richmond Ter Staten Island NY 10301


\(\color{skyblue}{\textrm{- Geolocating}}\)

library(tidygeocoder)

processed <- boro_halls_parsed %>% 
  geocode(address = pm.address, method = "osm", verbose = TRUE, unique_only = TRUE)

boro_halls <- boro_halls_parsed %>% 
  full_join(processed, by = c("pm.address" = "address" )) %>% 
  st_as_sf(coords = c("long", "lat"), crs = 4326)

boro_halls
hall street pm.address geometry
New York City Hall New York City Hall, New York, NY 10007 New York City Hall New York NY 10007 POINT (-74.00594 40.71274)
Bronx County Courthouse 851 Grand Concourse, Bronx, NY 10451 851 Grand Concourse Bronx NY 10451 POINT (-73.92332 40.82611)
Brooklyn Borough Hall 209 Joralemon St, Brooklyn, NY 11201 209 Joralemon St Brooklyn NY 11201 POINT (-73.9903 40.69276)
Queens Borough Hall 120-55 Queens Blvd, Queens, NY 11424 120-55 Queens Blvd Queens NY 11424 POINT (-73.8283 40.71378)
Staten Island Borough Hall 10 Richmond Terrace, Staten Island, NY 10301 10 Richmond Ter Staten Island NY 10301 POINT (-74.07609 40.64242)


\(\color{skyblue}{\textrm{- NYC Geocoding}}\)

#install.packages("remotes")
remotes::install_github("austensen/geoclient")
geoclient::geoclient_api_key("key")
geoclient::geo_search_data(boro_halls_parsed, street) %>% 
  transmute(address = input_location, 
            no_results, 
            bbl,
            bin = buildingIdentificationNumber,
            tract = censusTract2010,
            block = censusBlock2010,
            cd_name = communityDistrict,
            cd_num = communityDistrictNumber,
            puma = pumaCode,
            nta_name = ntaName,
            nta_num = nta,
            citycouncil = cityCouncilDistrict,
            assembly = assemblyDistrict,
            congressional = congressionalDistrict,
            latitude,
            longitude)
address no_results bbl bin tract block cd_name cd_num puma nta_name nta_num citycouncil assembly congressional latitude longitude
New York City Hall, New York, NY 10007 FALSE 1001220001 1079147 31 1021 101 01 03810 SoHo-TriBeCa-Civic Center-Little Italy MN24 01 66 10 40.71281 -74.00610
851 Grand Concourse, Bronx, NY 10451 FALSE 2024680001 2002869 63 3000 204 04 03708 West Concourse BX63 08 84 15 40.82627 -73.92297
209 Joralemon St, Brooklyn, NY 11201 FALSE 3001390001 3000256 11 1003 302 02 04004 DUMBO-Vinegar Hill-Downtown Brooklyn-Boerum Hill BK38 33 52 07 40.69248 -73.99049
120-55 Queens Blvd, Queens, NY 11424 FALSE 4022740002 4052812 216 1008 409 09 04111 Kew Gardens QN60 29 27 06 40.71332 -73.82858
10 Richmond Terrace, Staten Island, NY 10301 FALSE 5000070001 5000063 3 1004 501 01 03903 West New Brighton-New Brighton-St. George SI22 49 61 11 40.64244 -74.07533

\(\color{dodgerblue}{\textbf{Projection}}\)

\(\color{skyblue}{\textrm{- Checking}}\)

For vector data, use st_crs()

boro_halls %>% st_crs()
Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]


For raster data, just use crs()

canopy %>% crs()
CRS arguments:
 +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0
+ellps=GRS80 +units=m +no_defs 


\(\color{skyblue}{\textrm{- Reprojecting}}\)

Name Code
World Geodetic System 1984 (WGS 84) 4326
NAD83 / New York Long Island 2263

For vector data, use st_transform()

boro_halls_2263 <- boro_halls %>% 
  st_transform(2263)

boros_2263 <- boros %>% 
  st_transform(2263)


For raster data, use projectRaster()

canopy_2263 <- canopy %>% 
  projectRaster(crs = 2263)

\(\color{dodgerblue}{\textbf{Simplifying}}\)

For vector data, use st_simplify() (you need to transform out of lat/long first but can transform back after)

object.size(boros)
1372240 bytes
boros_simple <- boros %>% 
  st_transform(2263) %>% 
  st_simplify(preserveTopology = TRUE,
              dTolerance = 100) %>% 
  st_transform(4326)

object.size(boros_simple)
116320 bytes
str_c(round(100 * (116320 - 1372240) / 1372240, 2), "%")
[1] "-91.52%"
plot(boros$geometry)

plot(boros_simple$geometry)


For raster data, use aggregate()

manhattan_simple <- manhattan %>% 
  aggregate(fact = 5, fun = mean)
res(manhattan)
[1] 29.98979 30.00062
res(manhattan_simple)
[1] 149.9489 150.0031


ncell(manhattan)
[1] 619173
ncell(manhattan_simple)
[1] 24955


plotRGB(manhattan)

plotRGB(manhattan_simple)


\(\color{dodgerblue}{\textbf{Plotting}}\)

\(\color{skyblue}{\textrm{- Basic}}\)

ggplot() + 
  geom_sf(data = boros, aes(fill = population)) + 
  geom_sf(data = boro_halls, color = "red", size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) + 
  scale_fill_viridis_c()


\(\color{skyblue}{\textrm{- Basemap}}\)

library(ggspatial)

Options:


Bounding box to zoom:

boro_halls %>% 
  ggplot() + 
  annotation_map_tile(type = "cartolight", zoomin = 0) +
  geom_sf(color = "red", size = 1.5) + 
  coord_sf(xlim = c(-74.2591, -73.7002), 
           ylim = c(40.4774, 40.9162), 
           expand = FALSE) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank())


\(\color{dodgerblue}{\textbf{Analysis}}\)

\(\color{skyblue}{\textrm{- Centroids}}\)

st_centroid() creates a point at the center of each polygon

boros_centers <- boros_2263 %>% 
  st_centroid()

ggplot() + 
  geom_sf(data = boros_2263) +
  geom_sf(data = boros_centers, color = "blue", size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 


\(\color{skyblue}{\textrm{- Buffering}}\)

First check the projection to get the units you are working with (here, US survey foot, which is about 0.999998 feet)

boro_halls_2263 %>% st_crs()
Coordinate Reference System:
  User input: EPSG:2263 
  wkt:
PROJCRS["NAD83 / New York Long Island (ftUS)",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4269]],
    CONVERSION["SPCS83 New York Long Island zone (US Survey feet)",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",40.1666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",-74,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",41.0333333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",40.6666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",984250,
            LENGTHUNIT["US survey foot",0.304800609601219],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["US survey foot",0.304800609601219],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["US survey foot",0.304800609601219]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["US survey foot",0.304800609601219]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["United States (USA) - New York - counties of Bronx; Kings; Nassau; New York; Queens; Richmond; Suffolk."],
        BBOX[40.47,-74.26,41.3,-71.8]],
    ID["EPSG",2263]]

st_buffer() will buffer around a feature (though you cannot add units, so you have to know!)

buffers <- boro_halls_2263 %>% 
  st_buffer(dist = 5280)

ggplot() + 
  geom_sf(data = boros_2263) + 
  geom_sf(data = buffers, fill = "red", alpha = 0.5) + 
  geom_sf(data = boro_halls, color = "black", size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 


\(\color{skyblue}{\textrm{- Distances}}\)

st_distance() returns a matrix of the distance between each input to each output

boro_halls_2263 %>% 
  st_distance(boros_centers)
Units: [US_survey_foot]
         [,1]     [,2]      [,3]     [,4]      [,5]
[1,] 63944.00 25834.72  63119.49 29549.81 51982.247
[2,] 18452.27 21556.31 109797.41 66425.85 52004.377
[3,] 67591.47 31438.44  60926.31 21084.47 47928.670
[4,] 51683.39 44898.17 102397.11 41619.33  3501.901
[5,] 96115.36 57654.35  31051.68 35598.83 75287.273

More processing is needed to find, for example, the closest point

hall_center <- boro_halls_2263 %>% 
  st_distance(boros_centers) %>%
  as_tibble()

names(hall_center) <- boros_centers$boro_name

hall_center %>% 
  mutate(hall = boro_halls_2263$hall) %>% 
  gather(-hall, key = boro_center, value = distance) %>% 
  mutate(distance = as.double(distance)) %>% 
  arrange(hall, distance) %>% 
  distinct(hall, .keep_all = TRUE) %>% 
  full_join(boro_halls_2263, by = "hall")
hall boro_center distance street pm.address geometry
Bronx County Courthouse Bronx 18452.274 851 Grand Concourse, Bronx, NY 10451 851 Grand Concourse Bronx NY 10451 POINT (1005472 240261.9)
Brooklyn Borough Hall Brooklyn 21084.472 209 Joralemon St, Brooklyn, NY 11201 209 Joralemon St Brooklyn NY 11201 POINT (986940.8 191666.1)
New York City Hall Manhattan 25834.723 New York City Hall, New York, NY 10007 New York City Hall New York NY 10007 POINT (982603.5 198946.8)
Queens Borough Hall Queens 3501.901 120-55 Queens Blvd, Queens, NY 11424 120-55 Queens Blvd Queens NY 11424 POINT (1031849 199373.6)
Staten Island Borough Hall Staten Island 31051.678 10 Richmond Terrace, Staten Island, NY 10301 10 Richmond Ter Staten Island NY 10301 POINT (963132.8 173336.4)


\(\color{skyblue}{\textrm{- Area}}\)

st_area() returns the area of each polygon (make sure to check the CRS for the units!)

boros_2263 %>% 
  st_area() %>% 
  as_tibble() %>% 
  transmute(boro = boros_2263$boro_name,
            area = as.double(value)) 
boro area
Bronx 1187195250
Manhattan 636575581
Staten Island 1623637770
Brooklyn 1934173752
Queens 3040206185


\(\color{skyblue}{\textrm{- Length}}\)

st_length() returns the length of a line

library(tigris)

roads <- tigris::roads(state = "NY", county = c("Bronx", "Kings", "New York", "Queens", "Richmond"), progress_bar = FALSE) %>% 
  st_transform(2263)

bedford <- roads %>% 
  filter(str_detect(FULLNAME, "Bedford Ave"),
         !(LINEARID %in% c("110425921352", "110444047995"))) %>% 
  st_union()

ggplot() + 
  geom_sf(data = boros_2263) + 
  geom_sf(data = bedford, color = "yellow", size = 1) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank())

bedford %>% 
  st_length()
53673.32 [US_survey_foot]


\(\color{skyblue}{\textrm{- Dissolving}}\)

st_union() takes multiple entities and creates a single entity

mn_buffer <- boro_halls_2263 %>% 
  filter(hall == "New York City Hall") %>% 
  st_buffer(dist = 5000) %>% 
  st_transform(4326)

bk_buffer <- boro_halls_2263 %>% 
  filter(hall == "Brooklyn Borough Hall") %>% 
  st_buffer(dist = 5000) %>% 
  st_transform(4326)

ggplot() + 
  annotation_map_tile(type = "cartolight", zoomin = 0) +
  geom_sf(data = mn_buffer, fill = "red", color = "black", alpha = 0.5, ) + 
  geom_sf(data = bk_buffer, fill = "blue", color = "black", alpha = 0.5) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 

joined_buffer <- mn_buffer %>% 
  st_union(bk_buffer)

ggplot() + 
  annotation_map_tile(type = "cartolight", zoomin = 0) +
  geom_sf(data = joined_buffer, fill = "purple", color = "black", alpha = 0.5, ) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 


\(\color{skyblue}{\textrm{- Spatial Joins}}\)

st_join() returns everything, there is also st_intersects() and st_contains for more control

rand_points <- tibble(point = c("A", "B", "C", "D"),
                      lat = c(40.8845, 40.7202, 40.5905, 40.6366),
                      long = c(-73.7477, -73.9846, -74.0572, -73.7812)) %>% 
  st_as_sf(coords  = c("long", "lat"), crs = 4326) %>% 
  st_transform(crs = 2263)

rand_points %>% 
  st_join(boros_2263)
point geometry boro_code boro_name shape_area shape_leng
A POINT (1054014 261625.6) NA NA NA NA
B POINT (988518.8 201665.2) 1 Manhattan 636575908.328 359804.890798
C POINT (968363.5 154416.7) NA NA NA NA
D POINT (1044977 171282.8) 4 Queens 3040205398.17 900270.512975
ggplot() + 
  geom_sf(data = boros_2263) + 
  geom_sf(data = rand_points, color = "black", size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 


bk <- boros_2263 %>% filter(boro_name == "Brooklyn")

bk_intersect <- bk %>% st_intersects(roads)
bk_intersect <- roads[bk_intersect[[1]],]

bk_contained <- bk %>% st_contains(roads)
bk_contained <- roads[bk_contained[[1]],]

ggplot() + 
  geom_sf(data = bk) + 
  geom_sf(data = bk_intersect, color = "blue", size = 1) + 
  geom_sf(data = bk_contained, color = "yellow", size = 1, alpha = 0.5) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 


\(\color{skyblue}{\textrm{- Clipping}}\)

st_intersection() clips one entity based on how much of it intersects with another entity

boro_clips <- buffers %>% 
  filter(hall == "Bronx County Courthouse") %>% 
  st_intersection(boros_2263)

ggplot() + 
  annotation_map_tile(type = "cartolight", zoomin = 0) +
  geom_sf(data = boro_clips, fill = "purple", color = "black", alpha = 0.5, ) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 


\(\color{skyblue}{\textrm{- Bounding Box}}\)

st_make_grid() creates the smallest box possible containing all elements

boro_box <- boro_halls %>% 
  st_make_grid(n = 1)

ggplot() + 
  geom_sf(data = boros) + 
  geom_sf(data = boro_box, fill = "orange", alpha = 0.5) + 
  geom_sf(data = boro_halls, color = "black", size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 


\(\color{skyblue}{\textrm{- Convex Hull}}\)

st_convex_hull() creates the smallest shape possible containing all elements

hull <- boro_halls %>% 
  st_union() %>% 
  st_convex_hull()

ggplot() + 
  geom_sf(data = boros) + 
  geom_sf(data = hull, fill = "orange", alpha = 0.5) + 
  geom_sf(data = boro_halls, color = "black", size = 2) + 
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        rect = element_blank()) 


\(\color{skyblue}{\textrm{- Masking Rasters}}\)

mask() and crop() together clip a raster

canopy_2263 %>% 
  raster::plot()

canopy_2263 %>% 
  mask(boros_2263) %>% 
  crop(boros_2263) %>% 
  raster::plot()


\(\color{skyblue}{\textrm{- Overlaying Rasters}}\)

overlay() with a function can manipulate a raster

f <- function(rast) {
  rast > 60 & rast < 100
}

canopy_2263 %>% 
  overlay(fun = f) %>%
  raster::plot()


\(\color{skyblue}{\textrm{- Extracting Rasters}}\)

extract() pulls aggregated raster information

mean_canopy <- canopy_2263 %>% 
  extract(boros_2263, fun = mean)

median_canopy <- canopy_2263 %>% 
  extract(boros_2263, fun = median)

max_canopy <- canopy_2263 %>% 
  extract(boros_2263, fun = max)

min_canopy <- canopy_2263 %>% 
  extract(boros_2263, fun = min)

boros %>% 
  mutate(mean_canopy = mean_canopy,
         median_canopy = median_canopy,
         min_canopy = min_canopy,
         max_canopy = max_canopy)
boro_code boro_name shape_area shape_leng geometry population mean_canopy median_canopy min_canopy max_canopy
2 Bronx 1187193513.84 463868.937681 MULTIPOLYGON (((-73.89681 4… 2405464 13.062493 4.7337441 0 84.98292
1 Manhattan 636575908.328 359804.890798 MULTIPOLYGON (((-74.01093 4… 2736074 7.252479 2.1085426 0 76.45837
5 Staten Island 1623635734.34 325929.794027 MULTIPOLYGON (((-74.05051 4… 1694251 21.885133 15.4100775 0 87.58646
3 Brooklyn 1934174208.45 728195.197649 MULTIPOLYGON (((-73.86327 4… 1472654 3.652571 0.9797716 0 65.93524
4 Queens 3040205398.17 900270.512975 MULTIPOLYGON (((-73.82645 4… 495747 6.911773 2.9477946 0 86.00687

\(\color{darkblue}{\textbf{Leaflet}}\)


Used for creating interactive maps

library(leaflet)
library(leafpop)
library(leaflet.mapboxgl) # remotes::install_github("rstudio/leaflet.mapboxgl")

\(\color{dodgerblue}{\textbf{Basic Map}}\)

leaflet() %>% 
  addTiles() %>%
  addMarkers(lng = -74.000060, lat = 40.730910)

\(\color{dodgerblue}{\textbf{Zoom}}\)

Far (zoom = 10)

leaflet() %>% 
  addTiles() %>%
  setView(lng = -74.000060, lat = 40.730910, zoom = 10) %>% 
  addMarkers(lng = -74.000060, lat = 40.730910)


Close (zoom = 100)

leaflet() %>% 
  addTiles() %>%
  setView(lng = -74.000060, lat = 40.730910, zoom = 100) %>% 
  addMarkers(lng = -74.000060, lat = 40.730910)

\(\color{dodgerblue}{\textbf{Markers}}\)

boro_hall_points <- boro_halls %>% 
  st_transform(4326) %>% 
  st_coordinates() 

boro_halls <- boro_halls %>%  
  st_set_geometry(NULL) %>% 
  bind_cols(tibble(long = boro_hall_points[,1],
                   lat = boro_hall_points[,2]))


\(\color{skyblue}{\textrm{- Pin}}\)

boro_halls %>% 
  leaflet() %>% 
  addTiles() %>%
  addMarkers(lng = ~long, lat = ~lat)


\(\color{skyblue}{\textrm{- Colored Circle}}\)

boro_halls %>% 
  leaflet() %>% 
  addTiles() %>%
  addCircleMarkers(radius = 6,
                   color = "red",
                   fillOpacity = 0.8,
                   stroke = FALSE)

\(\color{dodgerblue}{\textbf{Labels}}\)

ny <- paste(sep = "<br/>",
            "<b><a href='https://www1.nyc.gov/'>City Hall</a></b>",
                 "City Hall Park",
                 "New York, NY 10007")

bx <- paste(sep = "<br/>",
            "<b><a href='https://bronxboropres.nyc.gov/'>Bronx Borough Hall</a></b>",
                 "851 Grand Concourse",
                 "Bronx, NY 10451")

bk <- paste(sep = "<br/>",
            "<b><a href='https://www.brooklyn-usa.org/'>Brooklyn Borough Hall</a></b>",
                 "209 Joralemon St",
                 "Brooklyn, NY 11201")

qn <- paste(sep = "<br/>",
            "<b><a href='https://queensbp.org/'>Queens Borough Hall</a></b>",
                 "120-55 Queens Blvd",
                 "Queens, NY 11424")

si <- paste(sep = "<br/>",
            "<b><a href='https://www.statenislandusa.com/'>Staten Island Borough Hall</a></b>",
                 "10 Richmond Terrace",
                 "Staten Island, NY 10301")

boro_halls <- boro_halls %>% 
  mutate(label = c(ny, bx, bk, qn, si))


\(\color{skyblue}{\textrm{- Popups}}\)

boro_halls %>% 
  leaflet() %>% 
  addTiles() %>%
  addMarkers(lng = ~long, lat = ~lat,
             popup = ~label)


\(\color{skyblue}{\textrm{- Hover}}\)

labels <- sprintf(str_glue("{boro_halls$hall}<br>{boro_halls$street}")) %>% 
  lapply(htmltools::HTML)

boro_halls %>% 
  leaflet() %>% 
  addTiles() %>%
  addMarkers(lng = ~long, lat = ~lat,
             label = labels)


\(\color{skyblue}{\textrm{- Formatting}}\)

leaflet() %>% 
  addTiles() %>% 
  setView(-118.456554, 34.09, 13) %>%
  addMarkers(lng = -118.456554, lat = 34.105,
             label = "Default Label",
             labelOptions = labelOptions(noHide = T)) %>%
  addMarkers(lng = -118.456554, lat = 34.095,
             label = "Label w/o surrounding box",
             labelOptions = labelOptions(noHide = T, 
                                         textOnly = TRUE)) %>%
  addMarkers(lng = -118.456554, lat = 34.085,
             label = "label w/ textsize 15px",
             labelOptions = labelOptions(noHide = T, 
                                         textsize = "15px")) %>%
  addMarkers(lng = -118.456554, lat = 34.075,
             label = "Label w/ custom CSS style",
             labelOptions = labelOptions(noHide = T, 
                                         direction = "bottom",
                                         style = list("color" = "red",
                                                      "font-family" = "serif",
                                                      "font-style" = "italic",
                                                      "box-shadow" = "3px 3px rgba(0,0,0,0.25)",
                                                      "font-size" = "12px",
                                                      "border-color" = "rgba(0,0,0,0.5)")))

\(\color{dodgerblue}{\textbf{Polygons}}\)

#First getting data from the census to plot (see section on census data)
library(tidycensus)
population <- get_acs(geography = "county",
                      state = "NY",
                      variables = c(population = "B01003_001"),
                      year = 2019, 
                      survey = "acs5",
                      geometry = TRUE,
                      progress_bar = FALSE) %>% 
  mutate(pretty_estimate = scales::comma(estimate),
         legend_value = estimate / 10000) %>% 
  st_as_sf()
pal <- colorBin("viridis", domain = population$legend_value)

labels <- sprintf("<strong>%s</strong><br/>Population: %s",
                  population$NAME, population$pretty_estimate) %>% 
  lapply(htmltools::HTML)

population %>% 
  leaflet() %>% 
  addTiles() %>% 
  addPolygons(fillColor = ~pal(legend_value),
              fillOpacity = 0.7,
              color = "white",
              weight = 1,
              opacity = 1,
              smoothFactor = 0.5,
              highlight = highlightOptions(weight = 5,
                                           color = "white",
                                           fillOpacity = 0.7,
                                           bringToFront = TRUE),
              label = labels,
              labelOptions = labelOptions(style = list("font-weight" = "normal", 
                                                       padding = "3px 8px"),
                                          textsize = "15px",
                                          direction = "auto")) %>% 
  addLegend(pal = pal, 
            values = ~legend_value, 
            opacity = 0.7, 
            title = "Population (10K)",
            position = "bottomright")

Return to Main Page